function [theta, J_history, theta_history] = gradientDescent(X, y, theta, alpha, num_iterations)

    m = length(y);
    % J_history = zeros(num_iterations, 1);
    theta_history = zeros(length(theta), num_iterations);
    % theta_s = theta;

    for iter = 1:num_iterations
        % theta_history(:, iter) = theta_s;
        theta_history(:, iter) = theta;
        % theta(1) = theta(1) - alpha / m * sum(X * theta_s - y);
        % theta(2) = theta(2) - alpha / m * sum((X * theta_s - y) .* X(:, 2));
        theta = theta - alpha / m * X' * (X * theta - y);
        % theta_s = theta;
        J_history(iter) = computeCost(X, y, theta);
        delta = theta - theta_history(:, iter);

        if (delta' * delta < 0.0000001)
            fprintf('alpha = %.2f, for loop end when iter = %d\n', alpha, iter);
            break;
        end

    end

end
